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ABSTRACT 


An investigation of a method for finding an optimum com- 
pensator for use in an autopilot is described. Given a struc- 
ture with a number of real poles and zeros which is steady- 
state decoupled, the vole-zero locations are chosen to give 
a system response as close as possible to a desired response. 
The desired response can be any one or a combination of sys- 
tem outputs. The computer 1s used to optimize the pole-zero 


locations by means of a function minimization program. 
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I. <GNTRODUCTION 


ties alcopLlOet 15 used fOr aircraft control to hold some 
variable of the aircraft's motion constant, such as speed, 
altitude or glide slope angle, in the presence of distur- 
bances. These disturbances can originate outside the air- 
craft, e.g. random wind gusts, or within the system such as 
electrical noise in actuators and control svstem components. 

The problem becomes more complex when the aircraft 1s 
viewed as a multivariable system, i.e. more than one input 
and output. For example, a movement of the elevator control 
will change the pitch angle and the speed, and a change in 
thrust likewise. It would be desirable for precise control 
to have one input cause a change in a corresponding output 
mmeon no other. This condition would require that the 
transfer function matrix of the entire system be a diagonal 
matrix, which allows much less freedom in designing for sta- 
Deity . 

This conditioncan be improved by requiring that the 
non-change in other variables (decoupling) occur only after 
the initial transients have died out, that 1s, once steady 
state has been reached. 

Methods have been devised to obtain steady state de- 
coupling by linear state-variable feedback and by cascade 
compensation with unity feedback (Huang). 

In particular, for the cascade compensaticn case, which 


1s probably more practical than state variablefeedback, Huang 





femme that the rank cond@etion for the plant matrix is no 
longer necessary as it was in the state variable feedback 
case and that pure integrators in the compensators and poles 
at the origin in the plant transfer function were helpful 

in providing steady-state decoupling. 

However, when it comes to designing the remainder of 
the compensators for stability or a desired transient res- 
ponse, one is left with a tedious root-locus design process 
which to some extent requires considerable experience in 
this technique. The problem is aggravated by the fact that 
the root loci for each channel are interdependent and this 
complexity increases with the number of inputs and outputs. 
Even when completed, there 1s no guarantee that the design 
is an optimum one. 

Mone intent of this thesis 1s to developa method of ap- 
plying the power of the modern digital computer and the tech- 
Nigues of linear programming to the process of compensator 
design. Essentially the process consists cf simulating the 
system in the time domain with the poles and zeros of the 
compensators as unknowns. By having a function minimization 
program minimize the integral-squared-difference between the 
response of the simulated system and a desired response, the 
requixed poles and zeros are determined. 

Cantalapiedra used a similar method with the same func- 


tion minimization program to determine the optimum number of 





poles and zeros te be used in a model of a high order system. 
In hiiis thesis.the integral squared difference between the 
system and model response curves was minimized for various 
configurations of the model and this so-called cost function 
served as a measure of how well the low-order model approx- 
imated the high-order system. 

Lima, in his thesis, used the same function minimization 
program to obtain the gains for a compensated ship controller 
in a model of two ships in underway replenishment. [In this 
case a desired track was to be followed by one ship with 
reference to the replenishing ship in the presence of surge 


and sway forces affecting both ships. 





Tt. DESC Or THE SYSTEM 


In this thesis, the longitudinal motion of the C-8A 
Augmentor Wing "Buffalo" aircraft presently under test and 
development at the NASA Ames Research Center was modelled. 
There are two reasons for this choice. First, a STOL air- 
craft such as the augmentor wing Buffalo, in the landing 
or take-off phase, operates at low speeds and steep flight 
path angles and in such flight regimes a strong artificial 
stability must be applied just to make the aircraft flyable. 
Second, a model of this same aircraft was used by Huang in 
his investigation of decoupling referenced earlier; thus the 
plant equations were readily available as well as the com- 
pensator transfer functions which he determined by the root- 
locus method. 

Only the longitudinal motion of the aircraft was in- 
vestigated since the principles arrived at for this mode 
would be equally applicable to the lateral-directional mode 
in a more complicated model. The aircraft 1s initially con- 
Sidered as flying straight and level at a speed Vo, = 126.7 
ft/sec when it is subjected to a reference pitch angle step 
Maput Of -0.25 radians at time = 0.0. These conditions 
apply to all aspects of the investigation which follows. 
Ideally, if the system were perfectly decoupled, the pitch 
angle would change to -0.25 radians and the speed and angle 


of attack would remain at their trim conditions. However, 


10 





as will be noted, there is a short period of transients 
which follows the initial step input which die out leaving 
speed and angle of attack at their trim Sen Glntone Ene the 
Ween angle at -0.25 radians. 

Figure II-l shows a schematic diagram of the system model 
consisting of the plant dynamics and the actuator dynamics and 
cascade compensators; one for each channel. The three inputs 
to the plant are Oy , % thrust (trim = 100%); O¢ , flap an- 
gle deflection from trim in radians and on , elevator deflec- 
tion from trim in radians. The three outputs from the plant 
are uu, the speed variation from trim in ft/sec, a , the an- 
gle of attack in radians, and 98, the pitch angle in radians. 
The feedback of outputs to inputs is to a certain extent ar- 
bitrary, however by inspection one could feedback an output 
to the input which had the greatest effect on it. In any case, 
the scheme adopted here is that used by Huang in his thesis. 

Figure II-2 shows a schematic diagram of the system 
model as designed by Huang. The compensator poles and zeros 
which he determined by root-locus techniques, the actuator 
transfer functions and the plant equations are also shown. 
The compensator poles at the origin, necessary for decoupling, 
Should be noted. The response of “the model to a pitch angle 
Seem input of -0.25 radians is shown in Figures II-3, I1I-4 
en@ I1-5 for speed, angle of attack and pitch angle respec- 
tively. These responses are used as reference responses 


during the course of the investigation. It will be noted 
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that these responses are somewhat oscillatory and of large 
amplitude. One major intent of this investigation is to 
improve (i.e., reduce) the frequency, magnitude and duration 
of these transient oscillations. 

It is assumed for this investigation that the structure 
of the compensators is fixed; that is, the numbers of poles 
and zeros in each compensator will be as designed by Huang. 
This is not to say that this is an optimum number. Indeed 
a greater number of poles and zeros in each compensator may 
improve the response over what is possible with the present 
configuration. A lesser degree of complexity and an inter- 
est in determining how much could be done with the present 
Cemerguration dictated the choice. Similarly, the gains in 
each compensator are considered as fixed at the values shown 
in Figure II-2, though there is no reason why thev could not 


be made a.further three variables. 


de? 





[Ifii. TEGUNI@QUE FOR SOLUTION 


Peewee COST FUNCTION 

The point in question for the investigation is to deter- 
mine whether and at what cost a given 3X3 plant can be made 
to follow a desired output response by determining the poles 
and zeros of cascade compensators of a given configuration 
using the digital computer. When given a curve of the de- 
Sired output response, the most obvious measure of how 
closely the model has matched it, is the integral-squared- 
difference between the deSired output and the actual output 
between t = 0.0 and some final time t . This so-called 


1g 


cost function J has the form: 


Xf 
n 
J = 2 
2. "i Yai ~ Yai) 4t 
1=1] 
where Yqi = the desired 1th output at time t 
a = the’ actug®. ith ouePpwt.atstime  t 
uae = Welgqiimngeractonr sor Sl the output 


The w, can be chosen to penalize a more important out- 
put deviation over a less important one or, as in the case 
of the present investigation, to give equal weight to all 
outputs which vary due to differences in units. It is not 
necessary that the cost function be a function of all outputs. 
It could be a function of the one of immediate importance 


with subsequent checking of the other outputs to follow. 


1 





Both types of cost functions were used in this investigation. 


The choice of te, is an important one. If te PEO 


long, the cost of the calculations in computer time becomes 
prohibitive, especially if the time steps used in integration 


are small. On the other hand if tt 1s too short, the re- 


si 


Suits may indicate a good following of the desired response 


mor t<t but the system may diverge and go unstable for 


5 aa! 
t>t,- A compromise was reached for this investigation 
choosing twotimes; a t, = 25 seconds, the approximate 


is 


settling time of the original theta response (see Figure II-5), 
and t, = 12.5 seconds: eThe ealeulatwens forecumve matchang 
were done over the shorter period and the results checked by 
Simulating the system and piotting the responses over the 
longer period. 

Thus it will be seen that J 1s a direct measure of how 
well the sytem output response curve matches a desired output 
response curve over the simulation time chosen. If J reaches 
a minimum it can be assumed that, barring the chance of a 
local rather than a global minimum, the "best" combination 
of poles and zeros has been found. Alternately, even though 
a minimum has not been found, a visual comparison of the two 


responses may indicate that the poles and zeros determined 


fee) Close enough to an optimum set. 


B. THE FUNCTION MINIMIZATION SUBROUTINE 
1. Generel 


Phew Eniere On minimization Subroutine chosen for tnis 


bo 





investigation is one which has seen some use at the Naval 
Postgraduate School in a variety of disciplines. At the 
time the present work was proceeding it was undergoing the 
incorporation of some improvements and the author had the 
benefit of this more up-to-date version. 

The subroutine is called BOXPLX and was programmed 
by R. R. Hilleary of the Naval Postgraudate School Computer 
@emcer. It will find the minimum of any arbitrary function, 
linear or non-linear subject to explicit constraints on the 
Yearebables or implicit constraints on functions of the vari- 
ables. It will handle a maximum of 25 variables, independent 
and auxiliary. With the subroutine must be appended two 
user-generated function subroutines to calculate the cost 
function and check for implicit constraint violations. A 
brief description of BOXPLX is included at Appendix A to 
allow a better understanding of some of the processes which 
follow and to guide further workers in this area. 

2. Constraints and Start Points 

The variables in BOXPLX are allowed to move within 
a feasible region (n-dimensional space, where n = number of 
variables) defined by upper and lower bounds on their values. 
These values are given as input to BOXPLX for any given pro- 
blem. 

It was assumed that all poles and zeros would be in 
the left half plane for reasons of stability and thus an 


absolute lower bound of 0.0 was set, but it may be necessary 
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Or desireable to set some positive vaiue as a lower bound 
depending on the results of a given run. 

The choice of upper bounds is more difficult. If 
they are set too high the feasible region in which BOXPLX 
searches for a minimum becomes large with a consequent 
heavy usage of computer time in reaching a minimum. If 
the bounds are too low, the "best" answer for a given run 
may be right on the bounds, necessitating a relaxation and 
further runs to proceed toward the minimum cost function. 
Some judgment 1S necessary in choosing these upper and lower 
bounds. 

The choice of the starting values for the variables 
is no less critical in avoiding wastage of computer time. 

A preliminary root-locus investigation of the system may 
give some idea of starting values which will be close to the 
"best" values, but it was decided that the real power of 
the method described herein would be in its application to 
the case where one had no idea where the answers were. Thus, 
the criteria mentioned above were not brought to bear too 
heavily on the problem. 

3. The Function and Constraint Evaluation Sub-Functions 

The function evaluation sub-function, FE(X,XDATA), 
is the routine which actually calculates the cost-function. 
For every set (vertex) of poles and zeros generated by BOXPLX 
is routine simulates the system response comparing it at 


every time step with the desired response, subtracting the 


“zak 





difference, squaring it and integrating the result from 
t= 0.0 to t = t. : 

KE(X) checks to see if there are any implicit con- 
straint violations in the generated vertex. Implicit con- 
straints may be any arbitrary function of the variables the 
user may desire, e.g. the product of two poles must be less 
than a given number. In this investigation no implicit con- 
straints were used. 

4. The Desired Response (XDATA) 

The response curve against which the system, with its 
eUement set of poles and zeros,is compared may be read into 
the programme as a set of data cards defining the amplitude 
of any arbitrary response at each time step of the simulation 
interval. Thus a curve which has no deterministic equation 
may be used to define a desired response. 

If the equation of a desired response curve is known, 
this equation could be included in function FE and the ampli- 
tudes computed simultaneously with the actual svstem response. 

In this investigation, data cards of the desired re- 
sponse were used, as being the most general method of spe- 
cifying the response. 

a oubroltine Graph 

This subroutine is appended to the main programme to 
plot the desired and actual response on a single graph when 
exit occurs from BOXPLX with the "best" set of poles and 


zeros. In addition, it can be made to plot other responses 


Zo 





which are not being compared with a desired response, to 


Smeek On their quality. 


SeCcONVERSION OF SYSTEM EQUATIONS TO STATE VARIABLE FORM 
i “General 

The system had previously been modeled by Huang 
uSing DSL (Digital Simulation Language) for which a program 
deck was available. However, since DSL is a language on its 
own, it was not possible to use it in its existing form; in 
an interactive mode in a Fortran program with the function 
minimization sub-routine. Thus, it was decided to first 
convert the system equations to state variable form so that 
they could be integrated by a built-in integration routine. 

Since the poles and zeros of the compensators will 
be the unknown parameters and it is intended to use the same 
quantity and configuration of these parameters as in Huang's 
thesis, there are nine unknowns; three poles and six zeros. 
The configuration of the compensators with these unknowns 
1s shown in Figure III-l, where the p, are the poles and 
the 2 are the zeros. 

BE. ' tne Plant “equatrons: 

The plant equations in Figure II-2 as excerpted from 
Reference 4 are almost in state variable form already except 
for the & term in equation (3). By substituting equation 
(2) into equation (3) and collecting terms the foliowing 


state equations result: 
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Substitution of the force and moment coefficients 
from Reference 4 into these equations results in the fellow- 


ing plant state equations: 


an = -.0670u + 22.20 =3220 + -06706, -6.126, (5) 
fe—eeoO40l7u - .76230 + .9479q ~- 0013556, ~.15276, --07526, (S) 
fe—eeoods785u - 1.826a-1.452q —- 0017326, +.099536,-1.94516,, (7) 
9 = q (8) 
3. The Actuator Equations 
ay Thrust actWator: The thrust actuator, in the no- 


fmmnon Of Figure Ii-2, has Y, as input and — as output 
meen transfer function: 


eval 
A + 1.98s + 1.21 





In differential equation form this becomes: 


5 § 2 
po leoeae + 12s, = 12k y, 


If we let Oy = § 


thrust actuator become: 


ta the state equations for the 


a Sig (9) 


6 = =. 210. — 1398s 


- cue aa (10) 


td 


b. Flaps actuator: The flaps actuator has input 


Y>5 and b¢ aS output With teansfer function: 
0 
s + 1.0 

In differential equation form this becomes: 
O¢ oF S¢ Fado 

Thus the state equation for the flaps actuator is: 
5. = —5¢ + Y5 (lar) 

c. Elevator actuator: The elevator actuator has 


input Y3 and @utput De with transfer function: 


NE We 8 
ae + Gos + lide. oO 


The differential equation for this transfer func- 


tron 1S? 
8 £66 we = 12 =  121.0y 
= e é 3 
Letting On = 8 og the state equations are: 
oe = bod ley 
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4. The Compensator Equations: 
In the notation which follows, the ij subscripts 
are defined as follows: the 1 Subscript refers to the 

thrust, flaps and elevator compensators as 1,2,3 respectively. 
The j subscript is the index of the coefficients in each 
compensator. The method of transforming the transfer func- 
tions into state equations when the resulting differential 
equation results in derivatives of the output is as set out 
in reference 5. In the notation of that reference, given an 


a order differential equation: 


x) (by + kee eee oa ae ap tht ante) = 
pou!” (t) i ee ie) a Bt ee 
ces 
= n 
where x) (4) a ae 
at” 


A set of state variables for this equation can be de- 


fined as: 
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Cr 
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x(t)-B u(t) 
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x(t) = Xz ft) Bu tt) 
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Bo = bo 

p= Bacay? 

By = bj-a,8)~-a58, 

BO = 2a ee seh earache = ay 197 a8, 


With this choice of state variable, the state equations are: 


X(t) = x, (t) +B u(t) 
R(t) = X(t) +B u(t) 
eet) = x(t) +8 ult) 
X(t) = a,x, (t)-a,_ x5 (t)-a,_ 5x3 (t)....-. ~a x, (t) +8 u(t) 


a. Thrust compensator: The thrust compensator has 


input e Output Yy and the transfer function is: 


10 
eee 

1 1000 (s+z,) (s+z.) (st+z,) 

=i S(stp,) (stp.) 


Multiplyimg out and transforming to a differential 


equation we have: 


y+(p,+p,)¥,+P Poy L908) O) + (21425424) e) + (2, 254252442)23) 2) 


+ (212524) 2, 
og, > (“Ahe ~ PYPp 4135 = 0 
by 9 = 1000 , b,, = 1000(z,+z,+2,) , b)5=1000(z,2,+2523+4 73) 
bi 3 = 10002 zz 
13 123 
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Thus the state variables are: 


Mee apes SY FS, 037% 7 F12%1 


Where 
B p=by 972000» By 4 =by Fy 92426 8p =P a 27 8112117? 10712" 
Bee 13 12°11 12 
and the state equations are: 
Wy = wo +B ey) (14) 
Wa = w+ Boe 1 5 ) 
Bee ~94.2%972 139" 39 13° ) 


The required output Vues when these equations are solved 


will then be: y, = w, + Boe 


b. Flaps compensator: The flaps compensator has 
input en output y., and its transfer: funetion is: 
Y> -400(stz,) 
C5 S 


Multiplying out and transforming to a differential equation 


we have: 
= -400e.-4002 ,e, 
Let Ao} = 0 
boo = -400 , boy = -4002, 


Thus for this first order subsystem, the state variable is: 


Wa Yo ~— Po0%2 
B = => - = = 6 
Where "59 = P20 400 , Py = Boi-Po9821 
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and the state equation is: 
=, = Px e (1) 


The required output Y> on solution cf this equation 
will be: 


ee 4° © o0°2 


©. Glevator Compensator: 


The input to the elevator compensator is e3 


and its output is Y3 and it has transfer function: 
¥3 -20(s+2,) (S+z_) 
e4 S(stp,) 


Multiplying out and transforming to a differential equation 


we have: 
V oe = Se EE : 
3+p3Y, 20 G+ (Z.+2,)e,+Z2.2-e, 
Let Az, = Pz + a3g5 = 0 
b36 ——— 2 ba) = ~20(2,+2Z,) ? b35 = ~202,2¢ 
The state variables are: 
95°23 °30°3 i yee algae 


Where Baq=b =—70e Bp 


30 BP ..5— 95578 Spares 


317317 3923278 32=P 3978 31 231-93 9432 


and the state equations are 


w= = W (18) 


is) 
We a 


gt P3123 


Sach aia 3 (1?) 


On solution of these equations the required output Y3 is 
obtained from: 


Y3 = WotB9e, 
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Pome Heedbackswand Cost Function Equations 
The feedback loops shown in Figure II-2 are included 
by simply adding three more equations rather than substituting 
in the state equations. This allows for more ease and flex- 
ibility in computer programming of the system model. The 


three equations are: 


a) = oe Ge 

C5 = To7 O 
= —¢ 

ae (3 


The general form of the cost function J discussed in Sec- 


tion above as applied to this system is: 


fe 
J Se =u)* + (ca - \2 + (8 = (2 at 
- > met ref ° ref ¥ 


ieee Complete System State Equations: 
Changing the variables in equations 5 to 19 inclusive 
to a consistent set cf variables as shown in Table III-1 the 
following complete system equations result: 


-.06/x,+22.2x.-32.2x,-6.12x.,+.067x 


256 


si IL Z 4 ] is) 
Xo = —-004017x)-.7263x,+.9479x,-.001355x,-.1527x5-.0732x_ 
X, = -003755x,-1.826x,-1.452x,-.001732x%,+.09953x 5-1. 945x_, 
Xy = X, 
Xe = X_ 

X¢ = ~1.21x,-1.98x,+1.21x,) )+1210.0e,) 
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Xo = “Xg+X)37400.0e, 

Xo = Xo 

hq = -121.0xg-6.6xg+121. 0x, 4-2420.00, 
St) 117" 11° 

my *12°? 12% 

ee) 2711 1112" 13° 
13 = 21% 

fy = * 157" 31%3 

His = ~437%15%30%3 

oe 1 1 

e9 = 575 

mo 374 


These equations and the defining equations tes 
the coefficients in terms of the free parameters (the poles 
and zeros) are coded in Fortran and appear in the sample pro- 
gram appended. The equations appear both in the function e- 
valuation subfunction to BOXPLX for calculating J and at)@ 
subroutine GRAPH for plotting the results. Further computer 
symbols used in this programme are explained in Dabsze Tfi-1. 
The computational flow chart showing the inter-relationships 


in the total programme is shown in Figure III-2. 
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TABES Tit - 1 

Variable Symbols 
Variable State Variable Computer Symbol 
X (1) 
XDOT (1) 
X (2) 
XDOT (2) 
X (3) 
XDOT (3) 
X (4) 
XDOT (4) 
X (5) 
X (6) 


=e & 
x“ OM 


* me OO 
An & & WW NY NY fF FY 


Ste 


Ore Om OF @D QO GQ R Re 
*% 


O 
cr 
Q 


XDOT (6) 
X (7) 
XDOT (7) 
X (8) 

X (9) 

X (10) 
agli: ) 
2) 

X (13) 

X (14) 

X (15) 

X (16) 
P} PA (1) 
P5 PA(2) 
P3 PA (3) 
PA (4) 
PA (5) 
PA (6) 
PA (7) 
PA (8) 


E £ OH & OF MH OF 
0 © & h @ o ct 
O 
© 
oF 
Mere e OO YD 
WwW N © 


at aS So see ee oe Peer 
> 


ON UW & W NY F 


a £ = 3 & 
= 
WN 


a3 





TABLE III - l (cont,) 


Variable Computer Svmbol 

Z 6 PA (9) 

ey El 

e, B2 

e3 B3 

ry R1 

ro R2 

r3 R3 

Upper bounds BU (1) 

Lower bounds BL (I) 

Starting values XS (I) 

Desired system output XDAPA(I), MDAPAT (1) 
Intermediate functions of parameters AIJ, BIJ, BETAIJ 
Time T 

Time step DT 
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oe yy OE INPUL LeU 
Desired response Nutbetmat trials Upper and lower bounds 
of system Number of points on variables 
Integration Interval ' Startinep, values 





MAIN PROGRAMME 


a BOXPLA : -FUNCTION KE(X) 


menerates points in bounded Checks for viola- 
epion and finds minimum of ron of Limpert 
a memet ton in that region constraints 


ee yp - OK KE=0,NOK=KE=-1 


ee 
System 

DUaee 

Equations 





an : oy 
FUNCTION FE CXDATA) 
KE<O?7 Yves SiMIIceS system waka 
given point, compares 
Wiley desired: OUlepm , 
calculates J, 












EX LT 
Minimum J is found 

Mimit on trials exceeded 
Unable to find feasible 
point which is better 




















SUBROLTINE GRAPH (X,XDATA) 
Plots system response using 
best values and desired 

system response, 






ews 


FIGURE ZTI1I-2 Computational Flow Chart 
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D> INVESTIGATIVE APPROACH 

Since it was not known how well BOXPLX would perform in 
this task of finding the poles and zeros it was decided to 
do all the initial investigation with the known original 
response with the poles and zeros as found by Huang. When 
this phase was completed, it was intended to specify an ar- 
bitrary but reasonable output response and attempt a matching 
of system response to this arbitrary one. 


As was noted earlier the only input used for all phases 


was a step input Kr =-70.25 radians corresponding to a change 
in the pitch angle 6-=-0.25 radians. Thus there are two out- 
puts, u and a which are not commanded by a corresponding 


input; and the question arises as to which output to make 

J a function of. This would depend on the situation in 
the real world. If the pitch angle were important due to a 
requirement to stay on a glide slope,then J should be a 
function of 86. If speed were important due to the distance 
mewceucndown, then J should be a function of uu. Ideally 
J should be a function of all three outputs, but it may 
not BO easy to specify just what they should be. The ap- 
proach used in the initial phase was to let J be a function 
of an output not commanded by an input (i.e., speed), then 
foe 2. be a function of pitch angle and finally a function 
of all three outputs. It should be remembered that in this 
initial phase, all three responses, and the poles and zeros 


which produced them, are known. The only other question to 


Si: 





Beecesolved is,that of the starting values and upper and 
lower bounds for the variables. Until the capabilities 
of BOXPLX had been checked it was decided to offset the 
Starting values approximately 25% from the known values 
(the so-called standard start point) with arbitrary close 
upper and lower bounds. These are shown in Table III - 2 


together with the actual values. 


Pole l Pole 2 Pole 3 
Actual 4.0 10.0 LOO 
XS Bao. ade) Q 9.0 
BU 50 L220 12200 
BL 320 §.0 8.00 
ZEGeo 1 Zero 2 Zero 3 zero 4 zero 5 Zero 
Actual LO 0.5 O25 Oss25 0.3 0.4 
XS O275 O55 0.335 O55 0.25 0.45 
BU 125 O27 0.7 0.7 0.4 O25 
BL O25 0:3 0.3 On3 0.2 O25 


Standard Starting Values 
CApgoe Lit = 2 

If BOXPLX proved capable of locating the correct poles 
and zeros to a desired degree of precision, the next step 
was to have it try with the so-called arbitrary start point. 
For this case it will be assumed that the pole and zero lo- 
cations are unknown and the starting values will all be set 
at XS = 10.0, the upper bounds at BU = 20.0 and the lower 
bound at BL = 0.0. 

Assuming success in this initial phase it was planned to 
then specify an arbitrary pitch angle response using Whiteley's 


Standard form of response to a step input with zero position 


o/ 


Jae 6 * ae ae © a eee 





error and 10% maximum overshoot with an appropriate w 
(Ref. 7). The speed response would then be checked and if 
"satisfactory", success will have been achieved. If not, 

a speed response using an overdamped second order impulse 
response would be combined with the pitch angle response to 


Meme a combined cost function. 


a6 





V. RESULTS AND COMMENTS 


A. INITIAL PHASE 

The initial phase, it will be recalled, consisted of in- 
vestigating the capabilities of BOXPLX when the poles and 
zeros of the desired response were known, i.e., the original 
poles and zeros as determined by Huang. It was also during 
this phase that various methods were used to reduce the 
amount of computer time involved in a run and to get an 
idea of how large the value of the cost function could be 
and still have an acceptably close response match. 

Table IV - 1 shows the results of various runs made with 
J a function of speed (1.e., a non-commanded output) and 
starting with the standard start point. The value of the 
cost function varies from J = .00036 to J = .052 . fhe 
response curves for the smallest value are so superimposed 
as to be indistinguishable and those for the largest value 
meemonmewn In Figures IV - 1, IV - 2 and IV - 3 for u , and 
a and @ respectively. Although the largest cost function 
is two orders of magnitude larger than the smallest one, the 
response matching is quite acceptable. Thus it appears that 
it is not necessary for BOXPLX to reach aminimum of the cost 
function to give gcod results and recognition of this fact 
could cut the computertime dramatically. 

the big differences in Table IV - 1 are in the computer 
time involved with each run. It was during this initial phase 


feat it was deerded to halve the simulation time to 12.5 


3 
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Semmes for the calculations in subfunction FE(X) and the 
results show that the computer time is also halved. Another 
ploy attempted was to have BOXPLX set up each complex with 
only n vertices instead of 2n vertices and this again 
halved computer time. However, this last method was found 
to work only when the starting values of the variables were 
close to the minimum and thus would be used only in the 
latter stages of an investigation. 

The next step was to investigate the effect of using the 
Seemme@ard Start point with J still a function of speed but 
enlarging the feasible region by setting the upper bounds at 
mueeemana the lower bounds at 0.0. BOXPLX found a cost func- 
tion of 0.05995 in 600 trials which is about twice the num- 
Pemweom trials it took to find the last entry in Table IV - 1 
which has a cost function approximately the same. The res- 
ponses with poles and zeros found on this run are shown in 
memes IV - 4, IV - 5 and IV - 6,for speed angle of attack 
and pitch angle respectively. - Again the results are quite 
acceptable. Another interesting result of this run is that 
the pole and zero values differ, in some cases appreciably 
both from those obtained from the last run in Table IV - l 
and from the actual values, and yet the responses are quite 
Similar. The values obtained from this run are shown below 
in Table IV - 2. A further mention of this sensitivity of 


the results will be made later. 
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25 


20 
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10 


time-seconds 


Figure IV = 6 
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Pole 1 Pole 2 Pole 3 


Found 4.409 0.437 EO 701 3 
Actual 4.0 FOTO OnO 

Zero l Zero 2 Zero 3 Zero 4 Zero 5 ZEXLO 6 
Found 0.664 Oe 25 O07 0.447 Ons 7 0.443 
Metual 1.0 OS Ces. Ono 0.3 O.4 
Parameters found, standard start, wide bounds, J=f{u) = .05995 


Raale: ly i= —2 


Next, the cost function J was made a function of the 
pitch angle output, 8 (the commanded ouput) again using the 
standard start point and the larger feasible region. in 600 


> and 


mates, taking 17 minutes BOXPLX found a J=6.8 x 10 
the results of the simulation of the poles and zeros found 
feemsnoOwn 1n Figures IV - 7, IV - 8 and IV - 9 foru ,a@ and 

8 respectively. Although the cost function as a function 

of 68 cannot be compared directly with the one which is a 
function of speed because of the scale factors involved it 

can be considered to a first approximation to be approximately 
three orders of magnitude smaller. Thus, this J asa func- 
tion of @ would be roughly equivalent to one as a évacuuae 

of Mm ef .068. But the important point to note is that the 
response match is much better than when J was a function of 
u and in fact the values found for the poles and zeros cor- 


respond to the actual ones with a maximum of 12% with another 


peee and the rest less than 2% Phus it would appear that 
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better convergence is obtained when J is made a function 
of the commanded output. 

Next came the real test of the capabilities of BOXPLX. 
Here it was assumed that the actual poles and zeros of the 
response are unknown and BOXPLX would start from a so-called 
Smeeitrary Start point in a large feasible region. For this 
Mueese, aS mentioned earlier, the upper bounds were all set 
at 20.0, the lower bounds at 0.0 and the starting values at 
10.0, as might be assumed by a reasonably competent eEngineer, 
examining this model. The cost function was made a function 
of speed and then a function of the pitch angle, as before, 
to check on the finding in the preceeding paragraph. 

feet, wath J = £(u), BOXPLX found J = 0.99016 in 736 
trials with the responses shown in Figures IV - 10, IV - 1l 
and IV - 12 for u,a and 6 respectively and found the poles 


and zeros shown in Table IV - 3 below. 


Pole l Pole 2 Pole 3 


Found ae 328 S2411 12.289 
metual 4.0 10.0 10.0 
zero l Zero 2 zero 3 zero 4 zero 5 Zero 6 
Found 0.229 On 722 1.342 Eee 17 0.084 12-130 
Actual 1.0 025 O25 0.5 O23 0.4 
Parameters found, arbitrary start, J = f(u) = 0.99016 
Table IV - 3 
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It will be noted that even though these values differ 
markedly from the actual ones, the u and © responses are 
mee but, there is a slowly decreasing steady state error 
in the response due to the very small value of Zero 5 which 
has almost cancelled the pole at the origin in the pitch an- 
gle compensator. 

In the hope of improving these results, the programme 
was rerun with the values in Table IV - 3 as the starting 
point, with the same upper and lower bounds (the bounds could 
be changed to narrow the feasible region, but it was decided 
to see how well the convergence went without doing this). 
firme turther 555 trials or a total of 1291 BOXPLX came up 


with a J = .08477 and the values shown in Table IV - 4 below. 


Pole l Pole 2 Pole 3 


Found 8.628 5.341 eS 


Actual 4.0 10.0 10°. © 

Zero l Zero 2 Zero 3 zero 4 zero 5 Zero 6 
Found 0.3000 L539 0.497 OeBe7 0 0.000 02912 
Actual ao) On '.5 OFS on 3 0.4 


Parameters found, arbitrary start, J = f(u) = 0.0848 


Table IV - 4 


The responses of thesystemwith these poles and zeros 
Mees NOwn in Figures IV - 13, IV - 14, and IV ~- 15 foru , a 


and 6 respectively. 
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fees tO be meted that, although J is smaller by an or- 
der of magnitude and the u anda responses are well nigh a 
Mempect Match, there 1S now a steady state error in the 9 
response due to a cancellation of the pole at the origin in 
the pitch angle compensator. This situation is not likely 
to improve with further trials with J a function of speed. 
This is the first indication of a worse response with more 
trials which will arise again later. 

The same arbitrary start point was used on the next run 
Beene J a function of 6 (the commanded output). This 
time, the results were more encouraging. In 900 trials with 

we 4.1 x io (again not directly comparable with the value 
when J 1s a function of speed) BOXPLX came up with the 


values shown it Table IV - 5. 


Pole l Pole 2 Pole 3 


Found 52259 8.708 o See 
Actual 4.0 Oe0 10.0 

zero 1 Zero 2 Zero 3 Zero 4 zero 5 Zero G 
Found 0.467 0.249 1. Gie2 0.447 Oa2725 0.374 
Actual ie. 6 OS 8.5 OS O53 0.4 
Seeemeters found, arbitrary start, J = £f£( ) = 4.1 x 10° 


Table IV - 5 
The responses for u, 4 and 9 are plotted in Figures 
IV - 16, IV - 17 and IV - 18 respectively. The match is near 
perfect for © and 9 and quite acceptable for u. This 


eo @irorcea the earlier assertion that better results 
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are obtained when J 1s a function of the commanded output 
(at least for this system). 

One new wrinkle which appeared during this last run was 
the case where BOXPLX was not able to find a feasible ver- 
tex at a given trial and restarted to compute at the best 
vertex. However, during the course of computation after 
this restart, 1t was again not able to find a feasible ver- 
tex and the best vertex was still the one at which it had 
restarted. Thus exit from the subroutine occured. The re- 
medy in this case was to change the value of R _, the first 
random number of the sequence generated in developing the 
initial complex of vertices. It must be admitted that this 
is probably a matter of luck as it is with all random numbers, 
but the results of the run speak for themselves. 

The last test in this initial phase was to make J a 
Function of all three outputs as an example of this admit- 
wedly artificial situation, since it 1s not likely that, for 
a given problem, one would know what all the outputs should 
look like, much less specify them. 

The problem arising here is to put weighting factors in 
the cost function so that the response variations from the 
desired response may be penalized depending on the importance 
of one or more outputs over others. If all are considered 
equally important which would normally be the case if they 
were specified it remains only to balance the cost function 


component variations due to scale factors. This was fairly 


64 





@aey in this case because of the data from past runs. The 
menue of J for the arbitrary start point when J was a 


function of speed was 0.9989046 x ie? and the equivalent 


Mer patch angle was 0.7685555 x oe. Dita! €un with J 
a function of all three outputs without weighting factors 


produced J = 0.9990991 x ia? Subtracting the first two 


from the last resulted ina J for alpha = 0.1176444 x er. 
Thus as a rough approximation, the speed component of the cost 
function was made 0.001, the angle of attack component, 0.1 
mn@ the pitch angle component, 1.0. 

mith this cost function, the best vertex after 1120 trials 


with J = .00042 was as shown in Table IV - 6. 


Pole l Pole 2 Pole 3 


Found b 2339 Taao3 Ona 7 
Actual 4.0 LOR D0 ag @ 

zero l Zero 2 zero 3 Zero 4 zero 5 Zero 6 
Found ios D id Det 1. 0918 0.641] 0.240 0.524 
Actual ike Q ).5 0.5 0.5 O23 0.4 
Parameters found, arbitrary start, J = f(u,a,6) = 0.00042 


Table IV - 6 


The system responses for these poles and zeros are plot- 
made in Figures IV - 19, IV - 20 and IV - 21 for u, a@ and 6 
respectively. The response matching is almost perfect al- 
though the pole and zero values differ considerably from the 


actual values. 
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me MRBITRARY DESIGN PHASE 

At this point, it was considered that the intricacies 
of BOXPLX and the factors affecting desired results were 
understood well enough to form a strategy by which a multi- 
variable system of a given configuration could be designed 
to have any time response desired to any input. 

The strategy was formulated as follows: 

a. Specify the time response and make the cost func- 
tion a function of the output which is commanded by the in- 
pec. 

b. Let BOXPLX try to match that output response for 


600 to 900 trials and examine plots of all outputs. 


c. If other outputs are reasonable, success has been 
achieved. If not, specify another output which may be import- 
Smee and make J a function of both outputs. This other out- 


put specification should be a compromise between what the 
ey eten seems capable of achieving and a reasonable desired 
response. 
d. Continue iteration until all responses are satis- 
Factory. 
Consequently, for this systemit was decided to specify 


the pitch angle response to a step input of -0.25 for r. as 


3 
Shown in Figure IV - 22. This particular step response was 
obtained from a state variable simulation of Whiteley's 


Standard form for a step response with a maximum of 10% 


Overshoot on zero position error (Ref. 7) with an Sa 8.75% 
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The value of Ww was chosen from an examination of the nat- 
ural time constant of tne system as revealed by previous data, 
which showed a rise time of approximately 0.5 seconds. [Ini- 
tials runs with a = 2.5 which gave a rise time of 1.75 se- 
conds showed that BOXPLX was having a difficult time trying 
to design the compensator for such a slow response. The 
values for Pole 3 kept ending up on the upper bounds, and 
they were raised twice to 30.0 and 40.0 before an ls B27 
was chosen. The computer program at the end of this paper 
shows how the response was simulated and cards punched to 
act as reference input to the function minimization program. 
This response is shown in Figure IV - 22. 

Thus a run was made with this response a reference and 
Meeemnetion of 8 . In 900 trials with J = .00156, BOXPLX 


found the poles and zeros shown in Table IV - 7. 


Pole l Pole 2 Pole 3 


3.669 14.644 19.890 

zero l HELO 2 zero 3 zero 4 zero 5 zero 6 

0.802 0239'7 0.034 2.054 Oe752 14220 
Meeameters found, arbitrary start, J = £(8) = .00156, ar- 


bitrary response 
Table IV - 7 


The response curves of the system with the above poles 
and zeros are plotted in Figures IV - 23, IV - 24 and IV - 25 
for u,4 and 9 respectively. The theta response matches 


the reference arbitrary response fairly well and the alpha 


as 





speed-ft/sec 


Jute 


14 


RZ 


10 





0 5 10 15 20 20 


time-seconds 


Figuee IV - 23 


TZ 





Pa 
| 
| 
_ 





—__ 
— a Sa 


© Ww) © Te) © Ww (>) 

N — om © ° © = 
® e ° ry () e 

= >) © © © 4 
! 


suetTper —- xoeQRQ}e JO aTbue 


sue 15 


eon 20 


Za 


20 


15 


10 


time-seconds 


Egogume IY = 24 


73 





pitch angle - radians 





E Swe 


a2 


<a 


0 is) 10 Jig 20 ZO 


time-seconds 


PIgure ye 25 


74 








Beeponse looks reasonable but there is an extremely long 
Meeclimg time for the speed to return to zero. This was 

due to the extremely small value of Zero 3 in the speed com- 
/emeator which almost cancelled its pole at the origin. 

In the hopes that more trials might improve the results 
the programme was run again starting from the values in Ta- 
fe EV - 7. After a further 900 trials or 1800 total, with 
J = .00088, BOXPLEX found the values give in Table IV - 8 


below: 


Pole l Pole 2 Pole 3 


peo) 3 6.635 14.933 

zero l Ze@ZO 2 zero 3 zero 4 zero 5 zero 6 

0.884 0.243 0.000 0.000 0.809 0.810 
Parameters found, arbitrary start, J = £(8) = .00088, ar- 


bitrary response. 
Treble IV - 8 


The results of asystemsimulation with these values are 
shown in Figures IV - 26, IV - 27 and IV - 28 for speed, an- 
gle of attack and pitch angle respectively. It can be noted 
that the match to the arbitrary theta response is excellent 
and that the short period oscillations both in it and in the 
alpha response have been reduced in amplitude and duration. 
But, the speed response now has a constant steady state error 
due to the fact that decoupling has been lost. It will be 
recalled that decoupling of this multi-variable system re- 


quired a pole at the origin in each of the compensators, but, 
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since Zero 3 in the speed compensator went to the lower 
Bound of 0.0, it effectively cancelled out this pole. 
Consequently, it was decided to change the lower bounds 
pee all the zeros to an arbitrary value of 0.25 and rerun 
fhe programme from the arbitary start point. This time, 
BOXPLX found the following values in Table IV - 9 after 


900 trials with J = .00180. 


Pole l Pole 2 Pole 3 


8.994 1 SieS 19 14.610 

Zero l Zero 2 zero 3 zero 4 zero 5 Zero 6 

te | 0.2918 2. 382 0.763 0.942 0.702 
Parameters found, arbitrary start, J = £(8) = .0018, arbi- 
meery response, BL = .25 on all zeros. 


Table IV - 9 


The response of the system with these values for the 
poles and zeros are plotted in Figures IV - 29, IV - 30, and 
IV - 31 for u,@ and 8 respectively. The results show 
some secondary overshoot on the @ response, some increased 
oscillations in the a response but the speed response is 
beginning to look like what one would expect from the con- 
ditions. 

To check if further trials would improve the results the 
programme was rerun for a further 900 trials starting from 
the values in Table IV - 9. The results of these 1800 trials 


mer J = .00099 are shown in Table IV - 10 below. 
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Pole l Pole 2 Fole 3 


7.007 8.231 m3 7 2 

Zero 1 Zeeo 2 zero 3 zero 4 zero 5 zero 6 

0.240 0.253 OuW567 0.240 0.842 0.828 
Parameters found arbitrary start, J = £(6) = .00099, arbi- 


@eery response. BL = 0.25 on all zeros. 


Tale IV - 10 

The responses of the system with these poles and zeros 
are shown in Figures IV ~ 32, IV ~ 33 and IV ~ 34 for uu, 

a and §@ respectively. AS waS mentioned before as being 
sometimes the case, these further trials did not improve the 
theta and alpha responses significantly but definitely wor- 
sened the speed response to one with a long settling time 
and some long period oscillation. 

Thus it was now time to backtrack to the previous run 
and specify a suitable speed response and make the cost func- 
tion a function of both u and6. The response chosen in 
concert with colleagues cognizant in this area was the one 
shown in Figure IV - 35. This response was obtained from a 
State variable simulation of the impulse response of a second 
Order overdamped system as shown in the computer programme at 
the end of this paper. 

After revising the function minimization programme to the 
final form shown at the end of this paper as a sample of all 
Such programmes, to include a balanced J asa function of 

u and 6 =, BOXPLX, in 1313 trials found the pole and zero 
values given in Table IV - 11. 
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Pole l Pole 2 Pole 3 


ie2). 236 222236 13.304 


zero l Zero 2 zero 3 zero 4 zero 5 zero 6 
1.444 @ 2310 2 ] 305 0.360 0.436 1. 210 
Parameters found, arbitrary start, J = £(u,8) = .04667, ar- 


bitrary response BL = 0.25 on all zeros. 


Table IV - ll 


ime results are shown plotted in Figure IV - 36, IV - 37 
mee iV - 38 for w,a and@ 6 respectively and the excellent 
results can be noted. This was also the best that BOXPLX 
could do since it exited after two restarts with no improve- 
ment in the best vertex. One need only glance at the ori- 


ginal response of the system to note the improvement obtained. 


Se oeNSLTIVITY OF RESULTS 

As was noted earlier, some excellent response matching 
has occured in the process of this investigation where, in 
the initial phase, the poles and zeros found by BOXPLX varied 
widely from those which produced the reference response. As 
a result there could be orders of magnitude difference in 
values of J and still the results could be acceptable. 
Thus, of itself, a "small" value for J is not necessarily 
indicative of a good response matching. 

A case in point is that where J was a function of all 
three outputs. The graphed responses were so superimposed 


as to be almost indistinguishable yet the value of Pole l 
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varied from its actual value by over 50% and that for Zero 3 
by over 100%. Both of these parameters were in the speed 
compensator and it might be said that these variations to- 
gether with those of the other three parameters in the same 
compensator might have interacted with each other to produce 
an overall acceptable response. Nonetheless, for all the 
parameters, the results show insensitivity for variations 


in the parameters of up to 25 and 30%. 
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¥V. CONCLUSIONS & RECOMMENDATIONS FOR FURTHER RESEARCH 
A. CONCLUSIONS 

The investigations carried out in this thesis result 
in the following conclusions: 

a. It is possible in a multivariable, steady-state 
decoupled system with cascade compensators and unity feed- 
back to design the compensators by computer methods to have 
the system match any reasonable specified response to a 
given input. 

b. Where computer time is a problem, there exist 
several methods, as outlined in this paper, for the function 
minimization subroutine used here to reduce the time taken 
to produce an acceptable response match. 

c. Where only a limited number of parameters are 
available, any arbitrary desired response should have time 
constants compatible with the natural time constants of the 
ayscem. 

d. Where possible, a cost function which is a func- 
tion of two or more outputs, when balanced, will produce 
better results than when it is a function of only one out- 
jour . 

e. Where the cost function is a function of only 
One output, the output to be chosen should be the one which 
is commanded by the input. 

f. The results show insensitivity of response to 


Variations in parameters of up to 25 and 303%. 
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B. RECOMMENDATIONS FOR FURTHER RESEARCH 
The following areas are recommended to extend the use 
of the design tool developed in this paper: 

a. The gains in the compensators, which were kept 
constant in this investigation, could be made free parameters 
to extend the flexibility of the design procedure. 

b. The numbers of poles and zeros in each compen- 
Sator were left at those developed in Huang's paper. An in- 
vestigation should be made to see if the response matching 
could be improved by increasing the numbers of poles and 
zeros, and by how many. 

c. The flexibility of the method would be further 
increased if provision were made to include complex poles 
and zeros as well as real ones. By uSing transfer functions 


2 
er the form Tu Ee ye where a and bb are real num- 


(S4p,) (S¥p,) 
bers, the parameters would be a combination of poles, zeros 
and coefficients. 

GQ. Military smecificationg for aircraft contain 
Specific values for limits on parameters of aircraft motion 
such as short period and phugoid frequency and damping. If 
the equations were altered to obtain these parameters of the 
aircraft response in terms of the compensator poles and zeros, 
a matching to military specifications could be obtained by 


insertion of a suitable implicit constraint evaluation sub- 


moutine to BOXPLX. 
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APPENDIX A 


DESCRLPTIOQN OF SUBROUTINE BOXPLX 


BOXPLX is a Fortran-63 subroutine which finds the mini- 
mum of a general non-linear function within a feasible region 
by the complex method of M.J. Box. Explicit constraints are 
defined as upper and lower bounds on the independent variables. 
iericit constraints may be arbitrary functions of the vari- 
ables. Two function subprograms to evaluate the cost function 
and the implicit constraints respectively msut be appended 
my the user. 

The complex method is an extension and adaptation of 
the simplex method of linear programming. Starting with any 
Ge feasible point (XS(I)) imn-dimensional space for n_ vari- 
ables, a "complex" of 2n vertices (sets of variables) is 
constructed by selecting random points (via a random number 
generator) within the feasible region. For this purpose, 

n coordinates are first randomly chosen within the space 
mBeunded by the explicit constraints. This defines a trial 
initial vertex, which is then checked for possible violation 
eeecne implicit constraints. If one or more are violated, 
the trial initial vertex is displaced half of its distance 
from the centroid ofpreviously selected initial vertices. 

This displacement process continues for up to 5n displace- 
ments until the vertex has become feasible. If not, exit 


@eeurs from the subroutine. 
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If an initial complex is established, the basic compu- 
tational loop is initiated. The current worst Vertex, 1.e., 
the one with the largest value of the cost function, is re- 
placed by its over-reflection through the centroid of all 
other vertices. (If the vertex to be replaced is considered 
to be a vector in n-Sspace, its over-reflection 1s opposite 
in direction, increased in length by a factor 1.3 and col- 
linear with the replaced vertex and the centroid of all other 
vertices). 

When this overreflection is not feasible or remains the 
worst vertex, it is displaced half-way toward the centroid. 
If NV displacements occur without a successful result, it 15s 
Over reflected through the best vertex, 1.e., the one with 
the smallest value of the cost function. If a successful 
result is still not achieved, the whole process restarts at 
the best vertex or the centroid, whichever has the smaller 
S@ec function. 

However, in most cases, all vertices converge to a point 
which is taken as a solution. Convergence is considered at- 
tained when 2n consecutive cost function values have been 
within 107° of a base value. 

The parameter list for BOXPLX is flexible and allows tne 
user to choose: 

a the number of variables 

b. the number of auxiliary variables (implicit con- 
Straints) 

©. te frequeney of printed output 


a. the number of trials allowed 
e. the first random number to be used 
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f. the starting values 
g. the upper and lower bounds 


If BOXPLX does not find a minimum cost function at the end 
of the number of trials specified it prints out the values 
for the best vertex reached up to that point. If further 
trials are considered necessary, these values should be used 
as the starting values for the next run. 

Obviously, the starting values selected must lie in the 
feasible region, i.e., between the upper and lower bounds, 
Or BOXPLX. will not even start. 

If calculation restarts at the best vertex after being 
unable to find a feasible vertex and the same situation oc- 
curs again, the "new" best vertex 1S compared with the pre- 
erous one. If it is better, calculations restart again at 
this "new" best vertex, otherwise exit from BOXPLX occurs 


with the previous best vertex being printed out. 
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